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Abstract. In this paper the two-dimensional structure of protoplanetary disks around Herbig Ae/Be stars is 
studied. This is done by constructing a self-consistent model based on 2-D radiative transfer coupled to the 
equation of vertical hydrostatics. As a simplifying assumption a grey opacity is used. It is found that the disk 
can adopt four different structures, dependent on the surface density distribution E(J?) as a function of radius, 
i.e. on radial- and vertical optical depth of the disk. For the case of high to intermediate vertical optical depth, the 
temperature and density structures are in agreement with the simple "disk with inner hole" model of Dullemond, 
Dominik & Natta (2001, henceforth DDNOl). At large radii the disk adopts a flaring shape as expected, and 
near the dust destruction radius (located at about 0.5AU for most Herbig Ae stars) the disk is superheated and 
puffs up. The region directly behind this "puffed-up inner dust wall" is shadowed, as predicted by DDNOl. For 
the case of intermediate to low vertical optical depth, but still high radial optical depth, the 2-D models show 
that the shadow can cover the entire disk. For such competely self-shadowcd disks the inner rim emission in the 
near infrared constitutes the dominant part of the SED, since the flaring component in the mid- and far infrared 
is suppressed by the self-shadowing effect. When the disk is optically thin even in radial direction, it becomes 
unshadowed again because the inner rim can no longer block the stellar light. Such disks have relatively weak 
infrared excess compared to the stellar flux. Finally, for disks that flare at intermediate radii, but become too 
opticaUy thin at large radii, the outer parts once again become shadowed. But this time the shadowing is caused 
by the flaring part of the disk, instead of the inner rim. The disk then consists of a bright inner rim, a shadow, a 
flaring part and finally a (dim) shadowed outer part. Different observational methods of determining the size of 
the disk (e.g. from the SED, from continuum mapping or from CO mapping) may yield different results. 
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1. Introduction 

Herbig Ae/Be stars are widely regarded as the interme- 
diate mass counterparts of the lower mass T Tauri stars 
(e.g. Strom et al. 



et al. 1997 
et al. 2001 
et al. 2003 



less s pheri cal geometry (di Fran cesco et al. 1994 ; Pezzuto 
Miroshnichenko et al. |l997|Jl999| ; Millan-Gabet 
and others favoring the disk-like picture (Vink 



Tuthill et al. 



Grinin & Rostopchina [1996| ; Natta et al. gOOl 
200ll ) 



1972; Palla & Stabler 1993; v.d. Ancker 



et al. 1998). But, while there is little doubt that the in- 
frared excess of classical T Tauri stars originates in a pro- 
toplanetary disk, the case for Herbig Ae/Be stars is not so 
clear. On large scales the disk-like nature of the circumstel- 
lar matter around Herbig Ae/Be stars is by now well estab- 
lished. Millimeter CO and continuum maps (Mannings & 
Sargent 1997), submm position- velocity maps (Qi pOOl| ) 
and imaging at vi sible and near- infrared wavelengths 
(Grady et al. 1999 ; Augereau 2001) clearly show rotat- 



Unfortunately, for some time disk models have failed 
to explain the SEDs of Herbig Ae/Be stars, and were con- 
sequently rejected by many. In particular the conspicuous 
bump around 3 microns in the SEDs of Herbig Ae/Be 
stars remained a mystery. Standard disk models such as 



those of Chiang & Goldreich (1997, henceforth CG97), and 
D'Alessio et al. (199^,1999) failed to explain this striking 
feature. Several explanations were suggested, ranging from 
FeO dust at 800 K (v.d. Ancker et al. |200(]| ) to accretion 
disks that are actively dis sipati ng only beyond a certain 
radius (Hillenbrand et al. 1992). None of these explana- 



ing flattened structures at scales of hundreds of AU. But 
on scales of tens of AU and smaller there seems to be con- 
flicting evidence, with some authors arguing for a more or 



tions were quite satis factory. Recently, it was recognized 
by Natta et al. (|200l| ) that this 3 micron bump may well 
originate from the inner rim of the dusty part of the disk. 
Since dust evaporates above about 1500 K, the inner parts 
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of the disk are free of dust. These gaseous inner parts have 
a much lower optical depth, and can even be entirely op- 
tically thin, dependent on the gas surface density. At the 
dust evaporation radius the dust forms a wall of 1500 K 
that is directly irradiated by the central star. This pro- 
duces an extra component in the spectrum that was not 
considered by the existing models of flaring disks. 

DuUemond, Dominik & Natta (2001, henceforth 
DDNOl) adapted the CG97 model to self-consistently in- 
clude this inner rim, and showed that the SEDs of Herbig 
Ae/Be stars can be naturall y exp lained in this way. In a re- 
cent paper (Dominik et al. 2002) th e sample of 14 Herbig 
Ae/Be stars of Meeus et al. ( 2001 ) was analyzed in the 
context of this model, and it was found that the SEDs of 
most stars were indeed consistent with the DDNOl pic- 
ture. 

According to this model, the inner rim of the flaring 
disk is much hotter than an ordinary flaring disk model 
would predict at that radius. This is because the inner rim 
is irradiated frontally rather than at a grazing angle. As a 
consequence, the inner rim is puffed up and casts a shadow 
over part of the flaring disk behind it. This shadow can 
extend from the inner rim, at about 0.5 AU, out to 5 AU 
or more. Outward of this shadowing radius the disk adopts 
the usual flaring shape as described by CG97. This part 
of the disk is responsible for the observed emission at long 
wavelengths. Dependent on the height of the inner rim, 
the shadow can reach so far out that the 10 fim silicate 
emission feature, produced by warm dust in the surface 
layers, is suppressed. This has been used by DDNOl as a 
possible explanation for the lack of 10 micron feature in 
several sources. 

Though succesful in explaining several features of the 
SEDs of Herbig Ae/Be stars, the DDNOl model was based 
on highly simplified equations. Among other things the 
structure of the inner rim and the shadowing of the disk 
behind it need closer theoretical examination. It is unclear 
from the DDNOl model what happens when the optical 
depth of the disk becomes too low to sustain flaring. Also, 
the DDNOl model was based on the assumption that, if 
the disk can flare outside the shadow, then it will. It is 
unclear whether perhaps in addition to these flaring disks 
also fully self-shadowed disk solutions exist. 

Because of the intrinsic 2-D nature of the problem, 
a closer theoretical study requires a full 2-D treatment 
of radiative transfer. This is done with a 2-D "Variable 
Eddington Tensor" solver. By coupling the radiative trans- 
fer to the equations of vertical hydrostatic equilibrium, 
the code solves the entire temperature and density struc- 
ture of the disk as a function of radius and vertical height 
above the midplane. As a simplifying assumption a grey 
opacity is adopted in this work. This is consistent with 
the disk consisting of large grains. The advantage of this 
simplification is that the results are more readily under- 
stood in terms of simple radiative transfer arguments. In 
a follow-up paper more realistic opacities and grain size 
distributions will be used, which will put us in a position 
to compare the results directly to observations. 



Using the 2-D disk structure code just described, the 
following issues will be addressed: 

1. Does the overall structure and SED of the disk agree 
with the much simpler model of DDNOl? Will there 
indeed be a shadowed region behind the inner rim, as 
predicted by DDNOl? How much will radial radiative 
diffusion heat this shadowed region, and how large is 
this region in radial direction? 

2. Under which conditions will the outer part of the disk 
collapse into the shadow of the inner rim, making the 
disk entirely self-shadowed? Can there be a bimodal 
set of solutions (self-shadowed and flared) for the same 
parameters? 

3. What is the structure of a fully self-shadowed disk? 
Will it be very cool and collapsed, or does radial ra- 
diative diffusion keep the disk still relatively warm? 

4. What will the outer (tenuous) part of a large disk look 
like: will it continue to be flaring, or will it sink into 
the shadow of the flaring part of the disk? 

The paper is organized as follows. In Sec. ^ the equa- 
tions are presented, and it is described how they are 
solved. In Sec. || a model of a canonnical flaring disk 
around a Herbig Ae star is described, and in Sec. ^ it 
is shown that self-shadowed disk can exist and what their 
structure looks like. In Sec. ^ it is investigated if multiple 
solutions exist for the same parameters. A description of 
the tenuous outer parts of a flaring disk is given in Sec. ^ 

2. The model equations 

The model is based on 2-D radiative transfer to calcu- 
late the dust (and gas) temperature in the disk. Once the 
temperature is found it is then used to integrate the equa- 
tions of hydrostatic equilibrium in order to find the density 
structure. The entire procedure is then iterated to obtain 
the full disk structure. Since the accretion rate in these 
disks is presumably very low, we do not include viscous 
dissipation in the equations. 

The model equations and the computational method 
used in this work are quite similar to the ones used 
for the 1-f 1-D (vertical) structure models presented in 
DuUemond, van Zadelhoff & Natta ( 2002 ). However, here 
this method is extended to full 2-D. Details of this method 
will be presented elsewhere (DuUemond in prep.). 

2.1. 2-D Radiative transfer 

Consider a gas -1- dust density distribution p{R, O) as a 
function of the spherical coordinates R (radius) and O 
(angle as measured from the polar axis). A grey opacity 

Kjy 1& adopted, and the scattering albedo is assumed 

to be zero. The star is assumed to be a point source. The 
flux of direct starlight F*{R, 9) at every point in the disk 
is then given by: 



F*{R,Q) 



L* 



AttR^ 



exp(-T(i?,e)), 



(1) 
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where t{R, Q) is the optical depth in the radial direction 
between the point {R, 6) and the central star, obeying the 
conditions: 

9r(i?,e) 



OR 



= p{R,e)K, 



(2) 



and r(0, Q) = 0. This direct stellar radiation is absorbed 
by the disk. The amount of energy per unit time and vol- 
ume that is absorbed in this way is: 



Q{R, e) = p{R, Q)kF*{R, e) . 



(3) 



By virtue of energy conservation this energy must then 
be reemitted as thermal (infrared) radiation by the dust 
in the disk. Dependent on the optical depth of the disk, 
this infrared radiation may be absorbed and re-emitted by 
the disk many times more, before it eventually escapes to 
infinity. It is this radiative diffusion process that allows the 
radiation to enter deep into the disk, and determine the 
temperature at every location, no matter at which optical 
depth. 

To solve this complex 2-D radiative transfer problem, 
the re-emitted radiation field is treated as a separate radi- 
ation field, which will be called the "reprocessed radiation 
field" from now on. The dust grains acquire a temperature 
such that they emit exactly the same energy per second 
as they absorb: 



TT An 



PkJ , 



(4) 



where J(i?, Q) is the frequency-integrated mean inten- 
sity of the reprocessed radiation field. The mean intensity 
J{R, Q) is defined as 



J = — [ i(n)dn. 



(5) 



where 1(^1) is the frequency-integrated intensity of the re- 
processed radiation in the direction f2. The first term in 
Eq.(^ accounts for the absorbed direct stellar radiation, 
while the second term takes into account the infrared ra- 
diation emitted by the disk itself. 

The radiative transfer equation for the reprocessed ra- 
diation is: 



as TT 



(6) 



which must be satisfied along every ray through the 
medium. Here s is the path length along the ray. The full 
2-D radiative transfer problem amounts to solving the cou- 
pled set of equations (||, ||, ^). The formal transfer equa- 
tion (Eq. ^ has to be integrated along all rays through the 
disk. In practice this means that a discrete (but large) set 
of rays is chosen which samples the reprocessed radiation 
field as accurately as possible. The integration of Eq.(^ 
along all these rays is done using the method o f "Ext ended 
Short Characteristics" (Dullemond & TuroUa pOOO| ) . 

Solving this set of equations is normally done using an 
iterative procedure from Eq.(||) to Eq.(^), to Eq.(||), and 



back to Eq.(^), until convergence is reached. This pro- 
cedure is known as "Lambda Iteration" and is the ba- 
sis of most current 2-D radiative transfer codes (includ- 
ing the original version of the program RADICAL described 
by Dullemond & TuroUa 2000). However, at large optical 



depth this method leads to convergence problems since 
information about the radiation field propagates only one 
mean free path per iteration. This can be fatal for prob- 
lems involving protoplanetary/protostellar disks. A rigor- 
ous solution to this problem is the "Variable Eddington 
Tensor" method, which couples the formal radiative trans- 
fer equation to the frequency-in tegra ted moment equa- 
tions (see e.g. Mihalas & Mihalas |l984|; Malbet & Bertout 
Stone, Mihalas fc Norman 1992; Dullemond, van 



1991 



Zadelhoff & Natta |2002D . RADICAL is now equipped with 
this method and can be used to solve 2-D radiative trans- 
fer problems at arbitrary optical depth with only a few 
iterations|^. 

2.2. Vertical pressure balance 

Once the temperature structure of the disk is determined, 
one can recompute the density structure by integrating 
the equations of vertical hydrostatic equilibrium: 



dP 

dz 



(7) 



where z is defined here as z = R{tt /2 — Q). Geometric cor- 
rection factors are neglected in the above equation, i.e. it 
is assumed that the disk is still thin enough that these 
factors are unneccessary. 

In order be able to solve Eq.(^, the surface density 
T,{R) = p{R, z)dz needs to be specified at every ra- 
dius. For a given value of S(i?), the vertical density struc- 
ture can be found by integrating Eq. (^ first for an initial 
guess of p{z = 0), and then renormalize it such that the 
required surface density is obtained. Having thus found 
the new density structure in this way, one can repeat the 
radiative transfer calculation for the next iteration step. 
Only when the new density structure is to within 10~^ of 
the previous density structure, this iteration procedure is 
terminated, and a solution is found. 

The function S(i?) is an input to the problem, and is 
entirely free to choose. This reveils an intrinsic weakness 
of models of passive disks: since one is not constrained by 
a global (constant) accretion rate, one has no theoretical 
constraint on this input function. In principle this intro- 
duces an infinite number of degrees of freedom. In practice 
this problem is less severe. Disks are presumably formed 
with a relatively smooth surface density profile. And the 
SED of a passive disk depends only weakly on small scale 
variations in the surface density. As will be shown below, 
significant changes occur only if large parts of the disk 
switch from optically thick to optically thin or vice versa. 



See web page http://www.mpa-garching.mpg.de/ 



PUBLICATIONS/DATA/radtrans/radical/. 
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2.3. Determining the SED 

Although the simphfying assumption of grey opacities is 
used in this paper, it is stiU useful to calculate the SED 
from the 2-D disk models. This is done by making images 
of the disk at a number of frequencies, and integrating 
these images to obtain the fluxes. The images are made 
using the ray-tracing capabilities of RADICAL. Care is taken 
that all spatial scales are sufflciently finely sampled. 

Computing the SED at a large number of inclination 
angles allows for a very powerful self-consistency check of 
the radiative transfer solution: the total outcoming lumi- 
nosity of the system should be equal to the luminosity of 
the star. At the end of each run RADICAL carries out such 
an energy conservation check. For all computed models, 
errors remained within 5% relative to the infrared lumi- 
nosity of the disk. 

3. A flaring disk around a Herbig Ae star 

As a first step, the structure of a "typical" flaring disk 
around a Herbig Ae star is computed. As mentioned above, 
the surface density function S(-R) is a parameter of the 
model. It is given as a power law of the type: 

= Eo (i?/AU)f (8) 

where p is the power law index. Furthermore an inner ra- 
dius Rin and an outer radius i?out is chosen, but it is veri- 
fied that the inner radius Rin produces a dust temperature 
that is consistent with the 1500 K dust evaporation tem- 
perature for silicates. For most models this is at 0.7AU, 
but for the optically thin model (model 3) a smaller inner 
radius is required. The parameters of the models are listed 
m Table |. The one discussed in this section is model 1. 

In Fig. 1^ the temperature and density structure of the 
disk is shown. The disk is considerably cooler at the mid- 
plane than at its surface. This is because the matter deep 
within the disk is shielded from direct stellar light, and is 
only heated indirectly by the reprocessed radiation field. 
The inner rim, on the other hand, is directly exposed to 
the stellar radiation field, and is therefore much hotter 
than the rest of the disk. Directly behind the inner rim the 
temperature drops strongly, down to values much smaller 
than the optically thin dust temperature at the same ra- 
dius. 

Between 6 AU and the outer edge the disk adopts a 
typical flaring shape. The flaring of the disk can be recog- 
nized by following the kink in the temperature gradient. 
The relative height (z/R) of this kink increases with ra- 
dius, showing that the surface height of the disk increases 
faster than R. This is the typical signature of flaring. 
This smooth temperature kink is the "surface layer" in 
which direct stellar radiation is absorbed and re-emitted 
as infrared radiation. Rather than a constant temperature 
(as in the simple model of Chiang & Goldreich and also 
DDNOl), this surface layer has a temperature gradient go- 
ing from the high temperatures of grains above the disk 




R [AU] 

Fig. 1. The 2-D temperature and density structure of a 
flaring disk: model 1 discussed in Sec.§. Contour lines rep- 
resent the temperature and are spaced 50 K apart. A small 
'-' symbol marks the 100 K contour. Grey scales are den- 
sity contours, logarithmically spaced such that 2 grey scale 
steps represent a factor of 10 in density. Dotted contours 
follow the grey scale contours, to aid the eye. The density 
contours are stopped at a gas density p = 10~^g/cm^ to 
avoid color crowding. The horizontal axis is the (spheri- 
cal) radius, in logarithmic scale. The vertical axis is the 
vertical height above the midplane measured in the angle 
7r/2 — 0, which is (to first approximation) equal to the 
vertical height z divided by the radius R. The equatorial 
plane is at the lower boundary of the figure. Note that 
in the coordinate system used here the stellar (primary) 
radiation moves along horizontal lines from left to right 
(i.e. at constant Q). Therefore, the shadow cast by the 
inner rim can be recognized as a steep jump in the tem- 
perature along an exactly horizontal line in the figure (at 
7r/2 — 6 = 0.18, ranging from the inner edge up to 6 AU). 



down to the lower temperatures of the grains deep within 
the disk. 

The structure of the disk between the inner edge and 
about 6 AU is different from the usual fiaring geome- 
try. In the figure one sees that the vertical temperature 
gradient suddenly acquires a very steep jump at about 
z/R ~ 7t/2 — Q — 0.18. This is a result of the shadow 
of the inner rim. Below this line there is no direct stellar 
radiation, while above it a dust grain would be exposed to 
the full stellar flux. In the outer regions (beyond 6 AU) the 
surface of the disk lies above this shadow line, and there- 
fore acquires the usual flaring disk shape. But at radii 
smaller than 6 AU the disk lies fully in the shadow, and is 
heated only by radial radiative diffusion. Note that, in the 
coordinate system used in Fig. stellar radiation moves 
along horizontal lines from left to right. 

The reason why the inner rim casts a shadow over the 
disk behind it is that the inner rim, being directly exposed 
to the stellar flux, is much hotter than the disk behind it. 
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Model 








rvert(50AU) 


'^rad 


P 




^out 


Type 


1 


2.0 


3.0 


10000 


1.9 X 10^ 


1 X 10^ 


-1 


0.7 


100 


Flaring 


2 


2.0 


3.0 


10000 


2.2 


5 X 10^ 


-2 


0.7 


100 


Self-shadowed 


3 


2.0 


3.0 


10000 


7 X 10"^ 


1 


-1 


0.3 


100 


Very tenuous disk 


4 


2.0 


3.0 


10000 


9.8 


4 X 10^ 


-1,-2 


0.7 


1000 


Flaring -I- shadowed 



Table 1. The parameters for the models in this paper. Column 2,3,4 are the stellar parameters: mass in units of Mq, 
radius in units of Rq and effective temperature in Kelvin. Column 5 is the vertical optical depth from = to = tt 
at a radius of 1 AU. Column 6 is the powerlaw index for the surface density such that t{R) = r(lAU) * (i?/AU)''. 
Column 7,8 are the inner and outer radius of the disk in AU. Colum 9 is for comments. 



0.40 



0.30 



0.20 



0.10 



0.00 




100 




Fig. 2. Same as Fig. |l|, but now computed using the stan- 
dard l-fl-D approach to radiative transfer. It is clear that 
the effects of the hot inner rim and the shadowed region 
are not present using this 1-l-l-D approach. 



and consequently has a higher pressure scale height. This 
hot inner rim is therefore puffed up. 

In Fig. ^ the same disk model was computed using the 
standard 1-l-l-D splitting, i.e. vertic al pla ne-parallel radia- 
tive transfer (see e. g. Ma lbet et al. 1991 ; Dullemond, van 
Zadelhoff & Natta 2002 ). This model nicely reproduces 
the flaring shape of the disk. But it fails to reproduce the 
heating of the inner rim and the shadowed region. This 
shows that, while 1-t-l-D models can be used reasonably 
well to model disks at radii larger than several AU, they 
cannot be used for the inner regions. 

Though the shadow is clearly seen in Fig. |^, in the 
shadowed region the midplane temperature does not drop 
very much. In fact, the temperature still increases as one 
goes towards smaller R. This shows that radial radiative 
diffusion manages to smear out radiative energy quite ef- 
ficiently, so that even matter in the shadow of the inner 
rim (all the matter with 7r/2 — 9 < 0.17 in this case) will 
still be relatively warm, though colder than what it could 
be if the surface of the disk were out of the shadow. The 
thermal emission from this shadowed region is somewhat 
suppressed, as can be seen in Fig. This suppression is 
less than what was expected by DDNOl. This is because 
DDNOl used a very approximate formula to estimate the 
effect of radial radiative diffusion. In the present 2-D cal- 



Fig. 3. The frequency-integrated intensity of the flaring 
disk (model 1) seen face-on, as a function of radius. This 
is computed by integrating the formal transfer integral 
(Eq. ^ vertically through the disk. The multiplied by R^ 
in order to weigh it with the emitting surface, so that it 
represents how strong the emission at each radius con- 
tributes to the SED. It is, so to say, the emission of the 
disk per constant interval in log(i?). 



culations the effect of this radial diffusion turns out to be 
stronger than predicted by DDNOl. 

From Fig. m it is also interesting to note that the emis- 
sion per unit of log(-R) is more or less constant for the 
flaring part of the disk (from about 3 AU onwards, up to 
shortly before the outer edge at 100 AU). Every radius of 
the flaring disk constributes equally to the SED. This is a 
feature that can be understood in terms of geometry. The 
surface height Hg of the disk goes roughly as Hs oc i?^/^ 
(see CG97). This means that the flaring surface of the disk 
captures stellar luminosity in a way that is proportional 
to i?^/^. This is only a slowly varying function of R. In 
Sees. ^ and ^ two cases will be discussed where the flaring 
geometry of the disk breaks down. And indeed, it will turn 
out that the "equal emission per log(i?)" will also cease to 
hold. Note that at the very outer edge of the flaring disk 
the emission drops somewhat (Fig. ^). This is a result of 
radiative leaking of energy through the outer edge, and is 
therefore a 2-D radiative transfer effect. 

In Fig. the SED of the 2-D model is shown at vari- 
ous inclination angles. One can see that for face-on incli- 
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Fig. 4. The spectral energy distribution (for grey opaci- 
ties) for model 1, computed at different inclination angles. 
The stellar atmosphere is included, but is approximated 
by a blackbody for simplicity. The inclination angles are 
measure from the pole (i.e. i = means face-on). Note 
that the vertical scale is linear. 



nations the 3-micron bump emission from the inner rim 
is very weak compared to the emission from the flaring 
part of the disk at longer wavelengths. This is because the 
inner rim is seen along the rim-surface, instead of perpen- 
dicular to it. At larger inclinations the inner rim emission 
gains in strength, while the flaring part fades. These are 
geometrical effects were predicted by DDNOl and seem to 
be confirmed here. It should be noted, however, that the 
strong supression of the near infrared bump at face-on in- 
clinations is a result of the perfectly vertical inner rim that 
is assumed in our model. In reality, hydrodynamic effects 
will presumably make this inner rim more round (as in 
Fig. |l^ below), which allows the near infrared emission to 
be seen also at face-on inclinations. As there is observa- 
tional evidence that even face-on disk show considerable 



near-infrared emission (e.g. Millan-Gabet et al. 1999), it is 
important to investigate such hydrodynamic "smoothing" 
of the inner rim in more detail in the near future. 

Images of how the disk actually looks like when viewed 
at a certain inclination angle arc shown in Fig. ^ in log- 
arithmic grey scale. A wavelength of 50 /im is chosen so 
that the entire disk radiates sufficiently brightly. In the 
upper-left image one sees the outer edge of the disk. It is 
evident that the surface layers radiate strongly, while the 
midplane layers are dimmer. As one zooms in, the intensity 
of the disk's surface gets stronger, and finally the bright 
inner rim appears. The inner rim is much brighter than 
the rest of the disk, as can be seen in linear grey scale (the 
lower- right image). The far side of the rim provides most 
of the flux. The inner rim therefore appears as an ellipse 
with one side brighter than the other. This asymmetry 
will presumably be less strong if one takes into account 
the smoothing of the inner rim by hydrodynamic effects, 
as discussed above. A strong shadowed region does not 




Fig. 5. Images of the flaring disk (model 1) at varying 
zoom factor. All images, with the exception of the bottom- 
right image, are in logarithmic grey scale. The bottom- 
right image is the same as the bottom-left image, but then 
in linear grey scale. From top-left to bottom-right the im- 
ages show the entire disk, zoomed in to the bright inner 
rim. The images are taken at 50 /Ltm. 

show up clearly in these images, dispite the fact that the 
shadow is indeed there (see Figs. |^ and||). Clearly the ra- 
dial radiative diffusion prevents the shadowed region from 
becoming too cold. 

4. Self-shadowed disks 

For disks that do not reach a certain minimal mass, it can 
happen that the entire disk will sink into the shadow of the 
inner rim. At no point the vertical optical depth is then 
large enough that the surface of the disk reaches out of 
the shadow. The disk (apart from the inner rim) is there- 
fore not directly heated by the stellar radiation, and will 
therefore cool down. Since the temperature determines the 
vertical height, the disk will therefore shrink even deeper 
into the shadow. The temperature and density structure 
of such a disk is shown in Fig.^ The shadow is now clearly 
seen to extend all the way towards the outer edge of the 
disk. Yet, the self-shadowed disk is not entirely cold, as 
can be more clearly seen in Fig. |^. Although no part of the 
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disk behind the inner rim can see the star directly, there 
is still sufficient indirect heating by the reprocessed radia- 
tion of the inner parts of the disk, and by radial radiative 
diffusion, that the self-shadowed disk remains relatively 
warm. 

The self shadowed disk has no hot surface layer: the 
temperature is more or less the same at every height above 
the midplane. Only at elevations significantly above the 
photosphere of the disk, where one gets in direct sight of 
the star, does the temperature rise to the optically thin 
dust temperature. But at those elevations the densities 
are already so low that one cannot speak of a hot surface 
layer anymore. 

In Fig. ^ the SED of the self-shadowed disk is shown, 
and compared to the SED of a flaring disk. Although this 
SED has been computed here only for grey opacities, it 
is interesting to note that the SED of the self-shadowed 
disk has a much steeper slope towards long wavelengths 
{lyF^ (X 1/^-^) than the flared disk (i^F^ cx which is 

a result of the steeper temperature slope of these disks. 
One may speculate whether such self-shadowed disks may 
be the explanation for the group II Herbig Ae/Be stars 
in the classification of Meeus et al.( ^OOl ). These sources 
show a weak far infrared emission, and have an overall 
SED shape very similar to the SED computed here for the 
self-shadowed disks. Moreover, the similar strength in the 
near-IR emission between self-shadowed and fiared disks 
is also observed between group I and group II sources. 
Yet, these group II sources have strong silicate emis- 
sion features, which might not be consistent with such 
self-shadowed disks. Also, mm observations indicate that 
group II sources are still quite massive, and it remains to 
be proven whether these disks are not too massive to be 
self-shadowed. To answer these questions, a similar disk 
model, but with more realistic dust opacities must be car- 
ried out. This will be one of the subjects of a follow-up 
paper (Dullemond & Dominik in prep.). 

5. Very tenuous disks 

When the optical depth measured along the equato- 
rial plane becomes smaller than unity, the height of the 
shadow drops to zero and the disk becomes unshadowed 
again. In this case not only the vertical optical depth Tvert 
is smaller than unity, but the optical depth along any ray 
through the disk is smaller than one. It can be consid- 
ered as a fully optically thin disk. Since such a disk has 
no surface layer which captures the stellar radiation (as 
in model 1), the disk is not a "flaring disk" as such. Yet, 
the density structure is still "flared" in another sense: the 
dimensionless pressure scale height hp = Hp/R increases 
outwards, as can be clearly seen from the density contours 
in Fig. 1^. This is easily understood, since the temperature 
of a grey dust grain goes as T cx 1/ \/R and the pressure 
scale height scales as Hp (x VtW, giving Hp/R cx R^^'^. 

Model 3 is a model with an equatorial optical depth 
(radial optical depth along the midplane) of Tjad = 1 . In 
Fig. p| the temperature and density structure of this disk 
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Fig. 6. The temperature and density structure of a self- 
shadowed disk: model 2 discussed in Sec. ^. Contour lines 
represent the temperature and are spaced 50 K apart. A 
small '-' symbol marks the 100 K contour. Grey scales are 
density contours, logarithmically spaced such that 2 grey 
scale steps represent a factor of 10 in density. 
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Fig. 7. The temperature at the equator as a function of i?, 
for the self-shadowed disk of model 2 (dot-dashed line) and 
for the more massive flaring disk of model 1 (dashed line). 
The temperature for a completely optically thin configura- 
tion is shown as a dotted line. Note that the temperature 
of models 1 and 2 near the inner edge (0.7 AU) exceeds 
the optically thin dust temperature. This is a backwarm- 
ing effect by the optically thick inner rim. 



is shown. It is clear that the temperature is only modestly 
influenced by the extinction at the equator, and that the 
overall density structure is that of a purely optically thin 
disk. The covering fraction of this disk is very small, and 
therefore the infrared excess of this disk is small compared 
to the steflar flux (Fig. |). 

The optically thin model presented here (model 3) may 
not be very realistic. In such very tenuous disks the tem- 
perature of the dust and the gas may start to decouple. 
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Fig. 8. The spectral energy distribution of the flaring disk 
of model 1 (dashed line), the self-shadowed disk of model 

2 (dot-dashed line), and the optically thin disk of model 

3 (dotted hne). 
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Fig. 9. The temperature and density structure of an al- 
most optically thin disk: model 3 discussed in Sec. ^. 
Contour lines represent the temperature and are spaced 
50 K apart. A small '-' symbol marks the 100 K contour. 
Grey scales are density contours, logarithmically spaced 
such that 2 grey scale steps represent a factor of 10 in 
density. 

Gas heating/cooling, dust-gas thermal energy exchange 
and dust drift should then be properly taken into account 
(see e.g. Kamp & van Zadelhoff 2001). 



6. Outer regions of flaring disks 

Generally it is assumed that the surface density S(i?) of 
a protoplanetary disk decreases as a function of radius. 
It can happen that a disk of the usual "rim -I- shadow + 
flaring" type will become too optically thin in the outer 
parts to maintain flaring beyond a certain radius i?turn- 
The Ty = 1 surface, where the stellar radiation is captured 
by the disk, in effect turns over and no longer maintains 



Fig. 10. The temperature structure of a large disk which 
becomes optically thin to stellar radiation at the outer 
parts: model 4 discussed in Sec. ^. Contours are 50 K 
apart. A small '-' symbol marks the 95 K contour. The 
dashed lines are also temperature contours, but with a 
spacing of 5 K, in order to make the structure at large 
radii more clear. Note that the radial range is larger than 
in the previous figures of this kind. 



the necessary flaring shape necessary to capture stellar 
radiation. The hot surface layer ceases to exist beyond 
that radius. The outer regions can now only be heated 
indirectly by the infrared radiation of the disk itself. This 
prevents the disk from collapsing to zero scale height. 

Model 4 is a model that has such a non-flaring outer 
part. To aggravate the situation so that the results are 
more clear, the surface density distribution is taken to be 
a broken powerlaw (p = — 1 inwards of 50 AU an p = — 2 
outwards of 50 AU). But the self-shadowed outer regions 
can occur also for single-powerlaw surface density slopes, 
though the dimming of the emissivity is then weaker. 

In Fig. ^ the temperature structure of this disk is 
shown. The shadowing effect is not very clear from the 
solid contours. That is why a number of finer-spaced con- 
tours are added in the figure (dashed fines). The shadow 
is apparent in the same way as it was for the shadow of 
the inner rim: a horizontally arranged temperature jump. 
But a better way to see the effect of the shadowing of 
the outer regions is to compute the emission of the disk 
per log(i?), as shown in Fig. |ll|. It is clearly seen that the 
continuum emission drops strongly in the shadowed outer 
regions, and will therefore only make a small contribution 
to the SED. 

This result may go some way towards explaining the 
differences that are often encountered when measuring 
the radius from SED fitting, mm continuum interfero- 
metric mapping an d CO line interferometric map ping 
(e.g. Mannings et al. |l997| ; |l997| ; Dominik et al. |2GG2| ). As 
the results of this model show, all of the radiation of the 
star is reprocessed within i?turn, since the disk is incapable 
of capturing radiation beyond that radius. So the SED will 
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Fig. 11. As Fig. |, but now for model 4. 



be dominated by the disk structure within i?turn, and the 
emission from the outer parts wiU be weak in comparison. 
Therefore the outer radius as inferred from SED fitting 
will yield i?out = -Rturn- But in interferometric maps one 
can see also the region beyond i?turn, even though it is 
much dimmer than the emission from the flaring parts of 
the disk. CO rotational lines, since they are density trac- 
ers, will be less senstive to this temperature drop. They 
may even be more easily detected in these outer regions, 
which may yield even larger effective disk sizes, unless CO 
is strongly depleted by freezing out onto dust grains. 

7. Bimodal solutions? 

In the paper of DDNOl the possibility was raised that 
two disk solutions might exist for the same disk param- 
eters, namely one where the disk has the rim-|-shadowed 
region -I- flaring part structure, discussed in Sec. |^, and one 
where, on the contrary, the whole disk is in the shadow of 
the inner rim, as discussed in Sec. |^. This is an important 
point because, if proven, it would mean that the actual 
shape of a disk depends not only on its properties such as 
mass and radius, but also on its previous evolution. 

The reasoning is as follows. Suppose we do not consider 
the process of radial radiative transport, and assume that 
the interior of the disk is only heated by radiation from 
the surface layer that diffuses downwards, but not side- 
ways to other radial shells. In such a scenario it is un- 
clear why the disk should always flare (and pop out of the 
shadow) at radii where it can. Suppose we could artifi- 
cally cool the flaring part of the disk down for a moment. 
Then the disk will shrink in height and retreat back into 
the shadow. Once in the shadow, it won't receive radia- 
tion from the star (the surface remains below the shadow 
line), and therefore the disk remains cool and stays in the 
shadow. Hence we would have found a second solution. In 
the DDNOl paper it was argued that presumably nature 
will pick the flaring solution only. Yet this remained to be 
proven. 



With the 2-D radiative transfer code one is now in a po- 
sition to do this, since 2-D transfer automatically includes 
the effect of radial radiative diffusion and indirect heating 
by radiation from the inner regions. The parameters of the 
disk are chosen such that it is only marginally flaring. The 
iteration of the disk structure computation was started 
from two different initial guesses for the vertical pressure 
scale height. One was a normal flaring disk {Hp cx R^^"^), 
and one was a fully self shadowed disk {Hp oc /Z~^/^). 

It was found that both initial guesses reach the same 
solution, which is a "rim -I- shadow -I- flaring" disk. This 
procedure was tried for a series of parameters and optical 
depths. In each case only a single solution was found. This 
means that the possibility of bimodal solutions seems to 
be excluded, at least within the present assumption of grey 
opacities. The conclusion is what one could call the "flar- 
ing disk principle" : a disk will flare whenever and wherever 
it can. Only if the disk cannot flare, because of too low 
optical depth, then it will remain self-shadowed. 

The question of bimodal solutions is also loosely linked 
to the question of stability. Based on a time-dependent 
version of the equations of CG97 for an irradiated flaring 
disk it can be shown that suc h disk s should be u nstable 
to self-shadowing (Dullemond ^000| ; Chiang ^OOOl) . Such 
disks would quickly develop ripples on the surface, which 
eventually cause most of the disk to collapse. However, 
these conclusions were based entirely on the highly sim- 
plified equations of CG97. It appears that 2-D radiative 
transfer effects will stabilize the disk (Dullemond in prep.). 
The radial thermal coupling between neighboring parts of 
the disk prevents the development of strong temperature 
gradients. 



8. Conclusions 

In this paper the structure of passive circumstellar disks 
was theoretically investigated. The problem was defined 
in a mathematically "clean" way: it is the problem of 
computing the temperature and density structure of ro- 
tating circumstellar matter around a star with a certain 
mass, radius and luminosity from basic principles of ra- 
diative transfer, radiative equilibrium and vertical hydro- 
static equilibrium. The disk parameters that went into 
the calculation were the inner and outer radius, and the 
surface density distribution as a function of radius. The 
only mathematical approximation made here was the re- 
duction of the hydrostatic equilibrium equations to 1-D 
vertical equations. From a physical point of view, many 
more approximations were made (related to dust-gas cou- 
pling, active accretion, dust opacities, etc). But these were 
necessary to keep the problem clear of uncertain physics 
for now. 

Four different kinds of solutions were found: a flaring 
disk, a self-shadowed disk, a transparent disk, and a flaring 
disk with self-shadowed outer region. These solutions are 
pictographically listed in Fig. ^ The numerical models 
described in this paper can be downloaded from a web- 
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C: Self-shadowed disk: D: Transparent disk: 




Fig. 12. Pictograms showing the four main kinds of solutions found. They represent a vertical cross-section of the 
disk, but are not to scale. A hashed area represents a shadow. In cases A, B and C the vertical optical depth may, 
under some conditions, drop below zero even though the radial optical depth remains much larger than unity. Case B 
is "zoomed out" to indicate that the second shadowing happens at large radii. In case D, the empty polygon is meant 
to show that even the radial optical depth along the midplane is small, meaning that the temperature everywhere is 
set by the optically thin dust temperature. The cases A, B, C and D are ordered according to decreasing optical depth 
(i.e. mass of the disk). It should be noted, however, that such an ordering does not always strictly apply since the 
powelaw index p of the surface density distribution also plays a role in determining the disk shape. 



site: www . mpa-garching . mpg . de/PUBLICATIDNS/DATA/ 
radtrans/grey2d/. 

The main conclusions are summarized as follows: 

1. A flaring disk around a Herbig Ae/Be star has a 
hot inner rim, a shadowed region behind it, and the 
usual flared geometry at large radii (Fig. |l2|-A) . These 
findings are in accordance with the predictions of 
DuUemond, Dominik & Natta ( |200l| ). But the effect of 
shadowing in suppressing the emission from the shad- 
owed region is not as strong as was predicted in that 
paper. 

2. Disks with intermediate to low vertical optical depth 
but high equatorial optical depth can become entirely 
self-shadowed (Fig. ^G). The SED falls off more 
steeply at long wavelength than for flaring disks. 

3. Disks with equatorial optical depths that are smaller 
than unity are un-shadowed again, since the inner rim 
can no longer stop the stellar radiation. These disks 
are fully optically thin (Fig. [T^D). 

4. The outer regions of flared disks can become shad- 
owed if beyond a certain radius the surface density 
becomes too low. This time it is the flared part of the 
disk that casts the shadow (Fig. [l^B). These outer 
parts do not contribute much to the SED, but may 
still be detectable using (sub-) millimeter interferome- 
ters. Measurements of the outer radius of a disk may 



therefore yield different results, depending on whether 
one uses SED-fitting, continuum mapping or CO map- 
ping. 

5. In the case of self-shadowed disks, and the shadowed 
outer parts of flaring disks, the usual 1+1-D approach 
to disk modeling breaks down, and so does the ap- 
proach used by CG97 and DDNOl. A full 2-D ap- 
proach, such as the one used in this paper, is then 
necessary. It might be that the SED is still relatively 
well described using a 1-l-l-D model or a DDNOl-type 
model up to the self-shadowing radius, but this re- 
mains to be proven using a 2-D model with more real- 
istic opacities. 

6. Bimodel solutions were not found. It seems that proto- 
planetary disks obey a kind of "flaring disk principle" : 
If the disk can flare, it will. Self-shadowed disks are 
therefore disks which cannot be made to flare. 

Many of these conclusions will presumably still hold 
when more realistic opacities are included. But this will 
be the topic of the second paper in this series. 
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